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An autoregressive model with a power-law type memory kernel is studied as a stochastic process 
that exhibits a self-afline-fractal-like behavior for a small time scale. We find numerically that the 
root-mean-square displacement A(m) for the time interval m increases with a power law as m“ with 
a < 1/2 for small m but saturates at sufficiently large m. The exponent a changes with the power 
exponent of the memory kernel. 


I. INTRODUCTION 

The Langevin equation is a typical stochastic differential equation used in statistical physics. Stochastic difference 
equations are also used to describe random processes that change stepwise by a unit of day, month, or year. The 
changes in daily maximum temperature and stock prices are examples of such random processes. The autoreCTessive 
(AR) model is a simple stochastic difference equation that has been used in statistics and signal processing Q. The 
AR model is expressed as 


k 

Xn — ^ ^ Xn—j T (1) 

1=1 

where Xn is a random variable, Oj’s (j = 1,2, ••• ,k) are parameters that represent the memory effect, and is 
assumed to be a Gaussian white random noise. In equilibrium statistical physics, the memory kernel and random 
noises are closely related by the fluctuation-dissipation theorem. However, such a relation is not usually assumed in 
the AR model, because thermal equilibrium is not considered. The AR models are simple linear models and used 
in various fields such as weather prediction and economic forecast. In the AR model, the time correlation function 
decays exponentially or is a sum of exponential functions. That is, the time correlation decays as C(m) = = 

Cj exp{ln(Aj)TO}, where cj’s are expansion coefficients and Aj’s are eigenvalues of the deterministic linear difference 
equation ii: 


— ajXn-j- 
j=l 


( 2 ) 


There are many linear difference models related to the AR model. The moving average (MA) model is expressed as 


i 

Xn — ^ + ^n- 

i=i 

A memory effect is included in the noised term. The autoregressive moving average (ARMA) model is a composite 
of the AR and MA models. The ARMA {k,l) model is expressed as 

k i 

Xn — ^ ^ kljXn—j T ^ ^ “f ^n- 
t=i i=i 


The difference in Xn can be expressed as 


Xn - Xn-l = (I - B)Xn 

using the lag operator B. The difference in rank d is defined as (1 — B^Xn- If two operators are defined as 

k i 

4>{B) = 1 - X 7r(H) = 1 + X 

i=i t=i 
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FIG. 1: Relationship between m and A(m) for the government approval rating in England. The dashed line is a line of 
A(m) oc with a — 0.375. 


the ARMA {k,l) model can be expressed as 


4'{B)Xn = TT{B)^n- 

The models using p, g of 0,1, 2 or 3 are used in most applications for fluctuating data. 

On the other hanchvarious types of fluctuating time sequence have been studied by many authors from the viewpoint 
of fractal physics [3,13. Recently, Honjo et al. have found a certain self-affinity in temporal fluctuations of the monthly 
government approval ratings of several countries [5|. Figure 1 shows A(to) = {{xn+m — as a function of to, 

where 0 < a;„ < 1 is the government approval rating in England between 1979 and 2013. (Original data are shown 
using percentage and is the value obtained by diving the data by 100.) A(to) changes as A(to) oc to“ with 
a = 0.375 for to < 25. However, A(to) tends to saturate for large to. In general, the self-afRne fractal with a < 1/2 
implies negative correlation. Many examples of self-affine fractals with positive correlation were reported in random 
interface problems il. In the normal diffusion, the root-mean-square displacement satisfies ((Ax)^}^'^^ ^ t°‘ with 
a = 1/2. The diffusion satisfying a yf 1/2 is called anomalous diffusion. The diffusion in random media sometimes 
exhibits anomalous diffusion. The exponent a yf 1/2 can be explicitly calculated for diffusions in some fractals M- 
Several mathematical models were constructed to generate self-affine fractals with a y^ 1/2. The fractional Brownian 
motion is a typical example Q. 

Some models related to the AR model can also describe random walks. The autoregressive integrated moving 
average (ARIMA) {k,d,l) model is defined using the difference operator B as 

</(R)(l - = 7r(R)C„. 

In the simplest random walk model, the time evolution equation is expressed as Xn = Xn-i + Cn- This is interpreted 
as an ARIMA model of p = 0, d = I, and q = 0. The operator {1 — B)‘^ for integer d can be expressed by the binomial 
expansion as 




d\ 


k=0 


k\{d-k)\ 


{-Bf 


If the rank number d is generalized to a non-integer, the model is called the autoregressive fractionally integrated 
moving average (ARFIMA) model [uHil. For the non-integer d, the binomial expansion is expressed using the 
gamma function as 


{1-Bf 


/ pcfc 

^F(fc+l)r(d-fc + l)^ ’ 


In the ARIFIMA model, the root-mean-square displacement A(to) obeys a power law A(to) ~ to“ with a y^ 1/2 for 
any to, and diverges for to ^ oo. 

In the problem of the government approval rating, A(to) needs to saturate for large to, because Xn changes between 
0 and 1. We would like to consider such a random process in which A (to) oc to“ for small m but saturates for large 
TO. Because {A(to)}^ = {{xn+m — XnY) = 2(x^) — 2{xn+mXn), the time correlation C{m) = (xn+mXn) decays as 
C{m) = C'(0)(1 — /3to^“). This decay law can be applied only for the time range where {A(to)}^ oc to^“ is satisfied. 
When TO becomes large, {A(to)}^ saturates toward a finite value and C{m) decays to 0. We do not yet understand the 
reason why government approval ratings exhibit such anomalous fluctuation, but the existence of obstinate persons 
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with long-term memory might be a reason for the anomalous fluctuation. Motivated by the anomalous fluctuations, 
we will study an abstract AR model with a long-term memory kernel in this paper. We study the simple AR model to 
reproduce the power-law behavior of A(rn) for small m. Our simple AR model is different from the ARFIMA model, 
in the sense that the difference operator B of the fractional rank d is not used, and A(m) saturates at large m. 


II. AUTOREGRESSIVE MODEL WITH LONG-TERM MEMORY 

As an AR model with long-term memory, we assume that a^’s decay in a power law. The model equation is written 
as 


Xn — (3) 

where A = j/Cq and Cq = j‘^)- Cg is the Riemann zeta function and it converges for <7 > I. That is, aj 

in Eq. (1) is set to be A/j'^ = j/{Cqj'^)- For 7 < I, the AR model Eq. (3) generates a stochastic process where 
the variance is finite. At 7 = 1, the variance (a;^) increases with n indefinitely like Brownian motions. This is 
because the deterministic linear difference equation without the last noise term. 





^n—j 


has an eigenvalue 1 at 7 = I, and Xn is not attracted to x = 0. 

Because the AR model is a linear equation, the solution can be obtained by the Fourier transform method. We 
consider a time series of length N. From Eq. (3), the Fourier transform x(wfe) = J2n=o Xne~'^‘^'‘^ for Wfc = 2TTk/N 
satisfies 


xicUk) = F{uJk)x{uJk) 

^N-1 

jj=0 


where F{u)k) = - x{uJk) is therefore solved as 


x{uJk) = 


The mean square displacement {A(to)}^ is evaluated as 


N/ 2-1 


1 - F[uk) ■ 


.2 n/2-1 

r A / \ T 9 4 V—> /I / \i9\r-. / xt 4(J V — > 1 — COS ( (jJh TTL ) 

{A(m)} =j^ ){l-cos(wfcm)} = — ^ \i-F(k)\^ ’ 


( 4 ) 


( 5 ) 


where denotes the strength of random noises. 

Figure 2(a) is an example of a random time sequence at g = 2, 7 = 0.98, and a = 0.5. The solid lines in Fig. 2(b) 
show the relationship between m and A(m) obtained by the direct numerical simulation of Eq. (3) for q = 2 and 
2.5 at 7 = 0.98 and cr = 0.5. The dashed lines denote the relationship between m and A(m) obtained by the 
numerical calculation of Eq. (5), which is consistent with direct numerical simulation. The difference between the 
data calculated using Eq. (5) and the direct numerical simulations is hardly visible in Fig. 2(b). A(m) obeys the 
power law A(to) oc for m < 100, and tends to saturate for large m. The alternate long and short dashed lines are 
the linear fitting in the double-logarithmic plot. The exponent a is evaluated as 0.223 for g = 2 and 0.325 for g = 2.5. 
The exponent a increases with g. 

The relation between the time correlation function C{m) and the parameters aj is known as the Yule-Walker 
equations [l^. If we assume a power law for the time correlation function, C{m) = C'(0)(1 — j3wP)^ there are three 
unknown parameters: C(0), j3^ and p. We can estimate the three unknown parameters from the time correlation C{m) 
at m = 0,1, and 2. Taking the average with respect to ^ after multiplying both sides of Eq. (3) by x„, <7(0) satishes 


C(0) 


{xD = 

m—1 


C{m) 


+ a 


.2 


( 6 ) 
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(a) (b) (c) 



FIG. 2: (a) Time sequence of Xn for q — 2, P = 0.98, and a = 0.5. (b) Solid lines denote the relationship between m and A(m) 
for the flnctuating time seqnences generated by Eq. (3) for q = 2 and 2.5 at /3 = 0.98 and cr = 0.5. Dashed lines denote A(m) 
calculated using Eq. (5). However, the difference between the solid and dashed lines is hardly visible owing to the overlap. Two 
alternate long and short dashed lines denote the lines of m“ with a = 0.223 and 0.325. (c) Relationship between q and the 
exponent a for 7 = 0.98. Rhombi denote the values evaluated from direct numerical simulations, and the solid curve denotes 
the result obtained from Eqs. (9) and (10). 


This is the first Yule-Walker equation for C{m) and a. Substitution of C{m) = C'(0)(1 — PmP) where p = 2a into 
Eq. (6) yields 


C(0) 


1 - 7 + iKq-pKq ' 


( 7 ) 


Similarly, by taking the average with respect to ^ after multiplying both sides of Eq. (3) by Xn-i and Xn- 2 , C'(l) and 
C(2) are expressed as 


C(l) = C(0)(l-/3)=C(0){7-dl/3^^^^^}, 

j=2 ^ 

°° ( i — 2]P 

C{2) = C(0)(l-/32^’)=C(0){-kl/3 + 7-kl/3V^:^-^}. 

3=3 ■' 


From the expressions of C'(l) and C{2), the following coupled equations are obtained, 

1-7 


P = 


i-7/C,E“2(j-i)Vi^’ 


l-7 = /3(2-f)-^/3E 


7 


(J - 2 )^ 




3=3 


( 8 ) 

(9) 

( 10 ) 


By the coupled equations Eqs. (9) and (10), p = 2a and P are solved at least numerically, and the variance (7(0) is 
calculated using Eq. (7). Rhombi in Fig. 2(c) show the numerically evaluated value of a and the solid curve denotes 
a = p/2 obtained from the solution of the coupled equations Eqs. (9) and (10) at 7 = 0.98. Fairly good agreement is 
observed. 

The parameter 7 determines the temporal range in which the power law is observed. Figure 3(a) shows the 
relationship between m and A(m) at 7 = 0.95,0.98, and 1 for q = 2.5 and cr = 0.5 obtained by direct numerical 
simulations of Eq. (3). Three lines overlap well for m < 20, and exhibit a power law of exponent a = 0.325. As 7 
decreases, A(m) saturates at small m. At 7 = 1, A(m) increases indefinitely. However, there is a cross-over in the 
power law of A(m), that is, A{m) ^ for small to, but A(to) ~ mP'^ for sufficiently large to. The dashed line 

in Fig. 3(a) is a line of wP'^. The exponent a = 1/2 implies that the statistical property of long-time fluctuations at 
q = 2.5 and 7 = 1 is equivalent to a simple random walk. We have evaluated the exponent a for sufficiently large to 
by changing g at 7 = 1. Figure 3(b) shows the exponent a as a function of q. The exponent a increases rapidly near 
q = 2 and becomes 1/2 for q > 2. 

In Fig. 4, we show the data of the government approval rating and A(to) for the AR model (3) at g = 3 ,7 = 0.95, 
and cr = 0.147. The characteristics of A(to) oc to“ with a ^ 0.375 for small to and the saturation for large to are 
reproduced by the AR model. 
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(a) 


(b) 



FIG. 3: (a) Relationship between m and A(m) at 7 = 0.95,0.98, and 1 for q = 2.5 and a = 0.5 obtained by direct numerical 
simulations of Eq. (3). (b) Relationship between q and the exponent a for 7 = 1 . 



FIG. 4: Relationship between m and A(m) for the government approval rating in England. The solid line is A(m) for the AR 
model at g = 3 ,7 = 0.95, and cr = 0.147. 


III. SUMMARY 

We have proposed an AR model with long-term memory obeying a power law. We have found that A(rn) increases 
with a power law of an exponent a < 1/2 and saturates for large m. In Sect. 1, we showed an example of a government 
approval rating that exhibits such anomalous fluctuations. We would like to investigate other stochastic processes 
in nature and social physics that exhibit a similar anomalous statistical property in the future. At 7 = 1 , A(m) 
increases m“ even for infinitely large m. The scaling behavior at large m is similar to that of the ARFIMA model. 
The construction methods of model equations are different; however, some mathematical relationship between the 
two models might exist, which is left for a future study. 
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